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ABSTRACT 

Current optical space telescopes rely upon silicon Charge Coupled Devices (CCDs) to 
detect and image the incoming photons. The performance of a CCD detector depends on 
its ability to transfer electrons through the silicon efficiently, so that the signal from every 
pixel may be read out through a single amplifier. This process of electron transfer is highly 
susceptible to the effects of solar proton damage (or non-ionizing radiation damage). This is 
because charged particles passing through the CCD displace silicon atoms, introducing energy 
levels into the semi-conductor bandgap which act as localized electron traps. The reduction 
in Charge Transfer Efficiency (CTE) leads to signal loss and image smearing. The European 
Space Agency's astrometric Gaia mission will make extensive use of CCDs to create the most 
complete and accurate stereoscopic map to date of the Milky Way. In the context of the Gaia 
mission CTE is referred to with the complementary quantity Charge Transfer Inefficiency 
(CTI = 1— CTE). CTI is an extremely important issue that threatens Gaia's performances: the 
CCDs are very large so that the electrons need to be transferred a long way; the focal plane 
is also very large and difficult to shield; the mission will operate at L2 where the direct solar 
protons are highly energetic (penetrating); and the science requirements on image quality are 
very stringent. In order to tackle this issue, in depth experimental studies and modelling efforts 
are being conducted to explore the possible consequences and to mitigate the anticipated 
effects of radiation damage. We present here a detailed Monte Carlo model which has been 
developed to simulate the operation of a damaged CCD at the pixel electrode level. This 
model implements a new approach to both the charge density distribution within a pixel and 
the charge capture and release probabilities, which allows the reproduction of CTI effects on 
a variety of measurements for a large signal level range in particular for signals of the order 
of a few electrons. 

Key words: instrumentation: detectors - space vehicles - astrometry - methods: numerical 
- methods: analytical - methods: data analysis 



1 INTRODUCTION 

We present a detailed physical Monte Carlo model of Charge Trans- 
fer Inefficiency (CTI) caused by displacement damage in irradi- 
ated CCD detectors. The development of the model took place in 
the highly challenging context of Gaia, a European Space Agency 
mission scheduled for launch in 2012. Gaia will operate for 5 
years in an orbit around the second Lagrange point (L2) (Perryman 
|et al.|200T||Lindegren et al.|2008| > and will measure the parallaxes, 
proper motions, radial velocities, and astrophysical parameters of 
over one billion stars. To do so, the satellite will constantly scan 
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the sky, observing the stars with two telescopes focussed on a sin- 
gle focal plane comprising 106 CCDs. The derived astrometric pa- 
rameters are highly sensitive to the precise image shape, and hence 
to the effects of CTI. 

Gaia will be subjected to the radiation environment at L2 
which is entirely dominated by protons emitted during solar flares. 
The energetic protons collide with and displace atoms in the CCD 
silicon lattice, leading to the creation of interstitial atom- vacancy 
pairs. The vacancies thus formed, combine, by diffusion, with other 
vacancies or impurities (e.g. oxygen, phosphorus, carbon atoms) 
present in the CCD as doping implants or due to pollution during 
the fabrication process. The impurity-vacancy complexes introduce 
energy levels in the semiconductor band gap that stochastically trap 
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and release the transferred signal carriers (electrons from the con- 
duction band in n-type devices). The time-dependent capture and 
release probabilities vary as a function of several factors; most im- 
portantly, the temperature, the local charge density distribution in 
the vicinity of the trap, and trap parameters such as energy level 
and capture cross-section. 

Based on the standard JPL model, the average proton dose re- 
ceived by the CCDs at the end of the 5 -year mission lifetime was 
originally predicted to be 4.14 x 10 9 protons cm -2 (with 90% con- 
fidence levels) for a launch in 2011 (Fusero 2007 ). Current space 
weather forcasts predict that the next solar maximum may be con- 
siderably less severe than average, so that the dose may be rather 
lower. However the sensitivity of Gaia to radiation damage is such 
that this in no way reduces the need to calibrate the effects. In addi- 
tion, the peak of the Solar activity is expected to occur in late 2013 
which means that Gaia will receive most of the total proton fluence 
early in the mission such that all data will be affected. 

Based on experimental studies led by the industrial partners in 
the Gaia project and independent analyses carried out within the 
Gaia science community, the CTI resulting from radiation-induced 
traps is expected to affect the mission performance by causing 
charge loss and image distortion. Mitigating those effects has been 
recognized as critical to achieving the mission requirements. Sev- 
eral aspects specific to Gaia contribute to the high impact of ra- 
diation damage on the mission. The large focal plane is difficult 
to shield and the CCDs will be exposed to most of the incoming 
particles. The required image location accuracy is extreme, e.g., 
the end of mission parallax error is required to be better than 25 
micro-arcseconds for a star of magnitude 15. The corresponding 
requirement on the residual image location error per CCD transit is 
0.01 pixels. However the image profile distortion induced by CTI 
has been measured to cause biases in the image location measure- 
ment of up to 0.17 pixels. Gaia will also study very faint objects 
down to magnitude 20 and at this signal level only a few electrons 
comprise the PSF core and the effects of trapping are poorly un- 
derstood. Gaia will scan the sky by continuously spinning around 
an axis perpendicular to the plane containing the telescope viewing 
directions (see Perr yman et al.|200"T||Lindegren et al.|2008| >. In or- 
der to follow the resulting motion of the stars across the focal plane 
and integrate the light during the transit, the CCDs will be oper- 
ated in Time-Delayed Integration mode (TDI mode). In this mode 
even fairly bright objects will remain faint for a part of their transit. 
Likewise, the sky background will form a gradient in the CCD par- 
allel direction and, due to the relatively short integration time, Gaia 
will not benefit from the potential trap filling effects of a bright sky 
background. 

These special operating conditions (high radiation dose, low 
signal level, low sky background and extremely high accuracy im- 
age location) demand a very high level of detail in the simulation 
of radiation damage effects, and preclude the use of models that 
assume instantaneous trapping within a certain volume that varies 
with the signal level (e.g., |Massey et al.|2010l|Rhodes et al.|2010| ). 
The Monte Carlo model presented here was developed to provide 
the required level of simulation detail. Our model thus simulates 
charge transfer at the electrode level and simulates the signal car- 
rier trapping thanks to a new approach to the representation of both 
the charge density distribution and the capture and release proba- 
bilities. These simulations are used within the Gaia project to study 
the effect of CTI on measurements, to generate simulated data with 
which to verify the future radiation damage mitigation algorithms 
and to obtain a better understanding of CTI itself. 



In the following sections we describe the relevant details of 
our Monte Carlo model and show that it can reproduce the experi- 
mental data obtained from irradiated CCDs operated in TDI mode. 



2 MODEL DESCRIPTION 

As described by Janesick ('2001), a CCD needs to perform four fun- 
damental tasks to generate an image: charge generation, charge col- 
lection, charge transfer, and charge measurement. The primary goal 
of our model is to simulate the effects of charge traps, induced by 
displacement damage in the CCD. These traps affect principally the 
third fundamental task of a CCD, the charge transfer, by stochas- 
tically capturing, and releasing charges during their transfer from 
one set of electrodes to another. The charge collection process can 
also be affected if a trap present in the CCD field free or depleted 
region captures a freshly photo-generated charge drifting towards 
the CCD buried channel. We chose not to take into account this sec- 
ondary aspect and focus on the signal charge transfer only. Thus our 
model simulates exclusively the transfer of charges present in the 
signal confinement region under each electrode in the image sec- 
tion and the serial register (Fig.[T^). The signal confinement region 
is simulated as a box (Fig. [I]}) of which the dimensions are defined 
by the manufacturing characteristics of the CCD, i.e. the width of 
the electrodes biased high in the transfer direction, the width of a 
pixel perpendicular to the transfer direction (serial direction) and in 
depth by the depletion of an electrode biased high with no electrons 
underneath. Fringing fields present at the edges of the signal con- 
finement region reduce the actual volume (Seabroke et al. 2008). 
To avoid any arbitrary assumption on the induced volume reduc- 
tion, the dimensions of the signal confinement volume remain set 
to the manufacturing characteristics, while, as we shall see later in 
Sections [23] and \%2A\ the actual distribution of the charge density 
within this volume is constrained by experimental measurements. 
This compensates to some extent our ignoring the fringing fields. 

Gaia CCDs ( [Short et al.|2 005) are custom made by e2v tech- 
nologies and referenced as CCD91-72. They are back-illuminated, 
full frame devices and incorporate a number of specific features 
such as a charge injection structure and 12 TDI gates to integrate 
bright stars over a shorter distance and avoid saturation. Each pixel 
also contains a supplementary buried channel (SBC), and an anti- 
blooming drain. In Table [TJimportant parameters of the Gaia astro- 
metric CCDs are summarized. Our model is capable of simulating 
the effect of charge injection and also takes the SBC into account. 
The TDI gates and anti-blooming drain are not explicitly modelled. 

2.1 Simulation process 

Figure [2]presents a simplified version of the whole simulation pro- 
cess. The first step consists of defining an input signal and speci- 
fying the CCD characteristics such as the number of pixels in the 
parallel and serial directions, the number of electrodes per pixel, 
the clocking scheme, the operating temperature and so on. 

Prior to the actual trapping and transfer simulation, empty bulk 
traps are randomly distributed across the CCD according to their 
specified concentration. Within each pixel, the traps are assigned 
a position in space by randomly generating coordinates within the 
signal confinement volume. If necessary the trap position can be 
kept fixed to repeat experiments with the exact same simulated 
CCD. The traps can belong to different species defined by the pa- 
rameters described in Table [2] according to the Shockley-Read-Hall 
(SRH) formalism. 



Gaia clocking scheme 
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a) Electron packet transfer 

Alternate clocking scheme 
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b) Pixel modelling 
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Figure 1. (a) depicts the transfer of several charge packets from pixel to pixel in a four phase device. A high voltage is applied on the black coloured electrodes 
(later referred to as 'biased high') while no voltage is applied on the empty ones (at rest). The ellipsoidal shapes represent the charge density distribution under 
the electrodes. The packets interact with the traps of different trap species (circle and triangle) during the dwell time At. The filled traps (filled black symbols) 
can release their charge (arrows) at any time. To illustrate the importance of an electrode level simulation, two different clocking schemes are shown; note the 
charge density distribution successively stretching and contracting in the alternate case compared to the Gaia case. 

(b) illustrates the signal confinement volume, a key point of our modelling approach, and puts it in the context of the CCD as well as the supplementary buried 
channel. The CCD pixel rows are referred to here as TDI lines. 



Parameter 


Value 


General 


Number of pixels (parallel x serial) 


4500 x 1966 


Number of light sensitive pixels 


4494 x 1966 


Pixel size (parallel x serial) 


10 x 30 /im 2 


Operational temperature 


163 ± 3 K 


Image section 


Number of phases 


4 


Transfer period 


982.8 ^s 


Pixel FWC 


190 000 e~ 


SBC FWC* 


~ 1300 e~ 


SBC size* (1 st CCD half) 


10 x 3 /im 2 


SBC size* (2 nd CCD half) 


10 x 4 /im 2 


Serial register 


Number of phases 


2 


Transfer period 


0.5 fis (avarage) 


Pixel FWC 


475 000 e~ 



Table 1. The e2v CCD91-72 parameters. SBC stands for supplementary 
buried channel and FWC for full well capacity. 

★ nominal values; the simulated and measured values differ significantly. 



Parameter 


Description 


E t (eV) 


Energy level in the semiconductor gap 


crt (m- 2 ) 


Capture cross section 


nt (traps/pixel or m -3 ) 


Concentration 



Table 2. The trap species parameters. 



The charge collection step corresponds to the generation of 
charges under the CCD electrodes. Photon detection from a light 
source can be described as a Poisson process, thus we create photo 
electrons in the CCD using a random generator with a Poisson dis- 
tribution and a mean equal to the expected number of collected 
photons within the integration time. As a consequence only integer 
numbers of electrons are generated so that the actual physical pro- 
cess is reliably simulated, including the variance in the number of 
electrons generated. This is important in the Gaia context since a 
considerable number of star images will contain only a few elec- 
trons per pixel due to the faint nature of the observed sources and 
the operating mode of the CCDs. In TDI mode, the light source mo- 
tion is synchronized with the CCD charge transfer rate so that the 
charge profile continues to build up as the image travels across the 
CCD. The transfer period is equal to the integration time at each 
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Figure 2. Top level diagram of the simulation process. 



step, which determines the number of photo-electrons generated. 
This implies that even for high signal levels the number of charges 
transferred will be very low, at least during the initial transfers. To 
simulate the background illumination (sky-light and scattered light) 
the corresponding background count rate is added to the expected 
number of collected photo-electrons before invoking the Poisson 
random number generator. Electronic charge injections (CIs) can 
also be simulated. In order to save simulation time any previous 
exposure of the CCD to a continuous level of background illumina- 
tion can be simulated analytically. This is done by pre-computing 
(cf. Section [24} the corresponding trap occupancy level of each of 
the trap species and filling the corresponding number of traps prior 
to the transfer of the signal of interest. 

The capture and release process is simulated as follows. From 
the trap parameters, the density of charges at the trap position (cf. 
Section [23] ), and the interaction time, it is possible to calculate the 
probability of capture or release (cf. Section [22| , for a specific trap 
relative to its state: empty or filled. If a trap is empty, the capture 
probability p is computed and a random number R generated; if 
R < p, then the capture is triggered and a charge is removed from 
the charge packet. Correspondingly, for a full trap, the release prob- 
ability is computed and a random number is generated; if a charge 
is released, it is added to the closest charge packet (cf. Fig. [T^). 
This procedure is repeated for each trap in the CCD, pixel column 
by pixel column. The interaction time (cf. Section [2^2} is defined 
by the amount of time a charge packet stays under a given set of 
electrodes. It defines the temporal resolution of the simulation and 
depends on the charge transfer period, the number of electrodes and 
the clocking scheme. CCDs with two, three, or four phases can be 
simulated, and any kind of clocking scheme applied. This facili- 
tates, for instance, testing the radiation hardness of different CCD 
configurations, taking into account the specific measurements to be 
carried out. 



During the charge transfer step, the CCD electrodes are biased 
high or set at rest according to the predefined clocking scheme. The 
charge packets are then redistributed under the next set of biased 
high electrodes. Trapping, transfer and charge collection (in TDI 
mode) are repeated until the last charge packet belonging to the 
input signal reaches the serial register, at which point the simulation 
ends. 

During read-out the signal charges are collected. It is also pos- 
sible to simulate the charge and release processes in the serial reg- 
ister by repeating the same procedure as for the image section and 
making sure the illumination history is respected, i.e. all pixels in 
a line are processed, ordered by distance from the output amplifier 
(Fig.Q>). 



2.2 Effective charge capture and release probabilities 

In the SRH formalism charge capture and release are described 
as decay processes. One can derive the charge capture and release 
probabilities as follows. First let us consider a number of filled traps 
bailed- in an infinitesimal time interval dt, the number of released 
charges is proportional to iVmied- The proportionality constant is 
the release rate r r . The number of filled traps as a function of time 
can then be derived: 



dN m] 



dt 

iVfiiled(t) 



— -r r iVmied : 
= ^Vfuii,oe~ rr 



(1) 



where iVf u ii,o is the number of filled traps at t = 0. The frac- 
tion of filled traps remaining after a time t is then statistically equiv- 
alent to the probability for any specific trap to remain filled after a 
time t: 



Nmedjt) 

iVfuii o 



: Pmied(£) = 1 -Pr(t) = e 



(2) 



where p r is the probability that a filled trap releases an electron 
within a time interval t: 



p r (i) = l-e- r ' 4 . 
The release rate constant is given by: 

1 v -Et/kT 

r r = — = XxcrtVthn c e . 



(3) 



(4) 



where r r is the release time constant, a t the capture cross- 
section, E t the trap energy level in the semiconductor forbidden 
gap, X the entropy factor, x me ne ld enhancement factor, T the 
CCD operating temperature, k the Boltzmann constant, and v t h the 
electron thermal velocity: 



3kT 

ml 



(5) 



with to* = 0.5 TO e the effective electron mass. The effective den- 
sity of states n c in the conduction band is: 



„ /2TTinlkT^ 3/2 



(6) 



where h is the Planck constant. 
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Likewise one can derive the probability that an empty trap cap- 
tures an electron within a time interval t: 



Pc(t) 



r c = — = crtVthrie 

T c 



(7) 



(8) 



where r c is the capture rate, r c the capture time constant, and n e 
the electron density at the trap location (cf. Section [23] ). 

The characteristic interaction time between the traps and a 
charge packet is the dwell time A t . It corresponds to the elapsed 
time between two charge redistributions, during which n e remains 
constant. The dwell time is proportional to the transfer period and 
depends on the selected clocking scheme. A t is greater in the CCD 
image area than in the serial register, since a complete TDI line 
(~ 2000 pixels) must be read out during a pixel to pixel transfer 
in the image area. For the Gaia CCDs, A t varies from several tens 
of micro- seconds (serial register) up to a forth of milli- second (im- 
age section). The release time constants for certain trap species can 
be as short as several hundreds of nano- seconds like the A centre 
(oxygen- vacancy complex). These 'fast' traps can lead to multiple 
capture and release events during a dwell time and play an impor- 
tant role in the CTI associated with the serial register. As a con- 
sequence we shall follow Lindegren (1998) and introduce effective 
probabilities of charge capture and release that take into account the 
possibility of multiple capture and release cycles within one time 
interval. We consider the probability pmied for a trap of unknown 
state to be filled after a time interval t. N is the total number of 
traps and A^mpty the number of empty traps: 



N = Nmi ed + N t 



empty 



(9) 

We follow the same derivation as in (eq.[TJ but now accounting for 
the empty traps 



dN m] 



dt 



r c Aempty ~ V r Afilled , 



1 dNmied _ dps 



N dt dt ' 

= TcPempty — ^filled , 

= r c (l — Pfilled) - ^filled • 



(10) 



In these equations p e m P ty is the probability for a trap of unknown 
state to be empty after the time interval t. One can then derive the 
following general solution, assuming that r c and r T are constants: 



Pfilled(t) 



+ Cexp[-(r r +r c )t] . (11) 



r r + r c 

It is now possible to derive the effective capture and release proba- 
bilities, p c and p T . If the trap is empty at t = then pmied(O) = 
and C = —r c / (r T + r c ), thus: 

p c = pmied(t) — ^ — (1 - exp [- (r r + r c ) t]) . (12) 

Tr + T c 

Similarly, if the trap is filled at t = then pmied(O) = 1 and 
C = +r c /(r r + r c ). Hence: 



l-p r = Pfilled(t) = 



r c + r r exp [- (r r + r c ) t] 
r r + r c 



r r + r, 



'— (l-exp[- (r T + r c )i\) 



(13) 



driven models. The volume driven models assume instantaneous 
trapping within a certain volume that varies with the signal level. 
A density driven model necessitates the computation of the capture 
and release probabilities for each trap — taking the charge density 
in the trap vicinity into account — regardless of its location (no 
trap is a priori ignored). The density driven model thus requires 
the evaluation of the charge density distribution as a function of 
the signal level and location within the pixel signal confinement 
region. The confinement region is defined by the electrodes biased 
high. As can be seen from the work ofjHardy et al.| f 1998) and more 
recently passey et al.| pOTO] ) and |Rhodes et al.f poTO| ), the vol- 
ume driven approach is fairly successful in explaining experimen- 
tal data, in particular HST data. However the Gaia operating condi- 
tions differ significantly from the HST ones. The Gaia CCDs will 
be operated in TDI mode. In this mode the exposure time equals 
the charge transfer period. Thus the charge-trap interaction time 
is significantly decreased and instantaneous trapping cannot be as- 
sumed anymore. Moreover, as already stated in the introduction, 
Gaia will deal with very low levels of background and source sig- 
nal. Trapping in these particular conditions was investigated for the 
first time in studies related to the Gaia mission. The density driven 
approach proved to be necessary to explain and reproduce experi- 
mental results ( |Short|2007l|Seabroke et al.|2008| ) which show that 
the CTI effects are modified by an extremely small level of back- 
ground light (of the order of a few photons per second). The volume 
occupied by these very few electrons should be negligible and thus 
prevent any trapping from occurring according to volume-driven 
models. However the Gaia experimental studies showed that very 
few electrons are capable of filling a significant amount of traps. 
This can be explained by a density driven approach. Indeed, as the 
background light constantly illuminates the CCD, the long effective 
charge-trap interaction time for these very few electrons compen- 
sates for the very small charge density in the vicinity of each trap 
and trapping then becomes likely to occur. Our simulation method 
is also particularly convenient for accurately simulating CTI effects 
over a wide range of signal levels and for CCDs in which extra dop- 
ing implants (such as a supplementary buried channel) modify the 
pixel potential and induce non-linearities between the CTI effects 
and the signal intensity (cf. Section |4|. 

The charge density distribution within the signal confinement 
volume, which varies strongly with the CCD architecture, is thus a 
key parameter in CTI modelling. Although the distribution cannot 
be directly measured, in principle it can be accurately determined 
for a specific CCD architecture and signal level by solving simul- 
taneously the Poisson equation and the charge continuity equation 
(e.g. |Seabroke et al.||2009| > in order to find a consistent electrode 
potential and charge carrier distribution. This requires a detailed 
knowledge of the CCD implant characteristics (i.e., the nature and 
concentration of the dopants). This information is often commer- 
cially sensitive and in order to keep our model flexible regarding its 
application to other cases than Gaia, we use an analytical descrip- 
tion of the charge density distribution which is roughly consistent 
with the modelling results by Seabroke et al. ( 2009). This analytical 
description consists of a normalized Gaussian function in the three 
space directions of which the complexity increases with the num- 
ber of the CCD potential characteristics included. The distribution 
parameters are listed in Table [3] The density is defined as: 



2.3 Charge density distribution modelling 



The CTI effects model we describe in this paper is a density driven 
model to be contrasted with the more commonly used volume 



n e (x, y,z) = S x p(x, y,z) = S x p(x) , (14) 
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Parameter 


Description 


S(e~) 


number of e — in the signal confinement volume 


n sa t (m -3 ) 


density saturation level 


Signal confinement volume 


£max (m) 


parallel width 


2/max (m) 


serial width 


2max (m) 


depth 


V (/im 3 ) 


signal confinement volume x ma x X y m ax X z ma * 


Buried channel regime 


(j x (m) 


parallel distribution standard width 


(Ty (m) 


serial distribution standard width 


(J Z (m) 


depth distribution standard width 


xo (m) 


distribution centre parallel coordinate 


yo (m) 


distribution centre serial coordinate 


z (m) 


distribution centre depth coordinate 


Ssat(e-) 


Signal level at which the buried channel saturates = FWC 


Supplementary buried channel regime 


Cy, SBC (m) 


parallel distribution standard width 


o-z, sbc (m) 


depth distribution standard width 


2/o, sbc (m) 


distribution centre serial coordinate 


*o, sbc (m) 


distribution centre depth coordinate 


#sbc (e~) 


Signal level at which the SBC saturates = FWCsbc 



Table 3. Charge density distribution parameters 



where 



with 



Hence 



exp [-|(x-x ) T C- 1 (x-x )j 
= — , ^3 |r , 1/2 1 , d5) 



(v^F) |C 



G% 

(Ty 

ol 



(16) 



exp 



p(x) = 



i(( £ e) 2 +( 3£ if) 2 +( £ e) 2 ) 



27r) a x a y a z 



(17) 



The supplementary buried channel (SBC) corresponds to an 
additional doping implant, which generates a deeper potential well 
and narrows the charge distribution so that electrons are transferred 
through a smaller volume of silicon and encounter fewer traps. In 
the Gaia CCDs, the SBC potential well collapses at signal levels 
of a few thousand electrons. As a consequence, the SBC only im- 
proves CTI in the small signal regime. To implement the SBC, 
we simply introduce a second, low signal regime with a different 
charge density distribution and ensure a smooth transition between 
them so that (eq.[l7} becomes: 



exp 



(^f) 2 +(^) 2 +(^y 



(18) 



where the parameters indicated with a ★ now vary as a function of 
signal level: 

P* = p(l- e~ ^fe) + P S Bce"^fe , (19) 
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Figure 3. Example of a charge density distribution as it is simulated by the 
presented analytical model: the charge density at the center of the charge 
cloud (in the parallel direction and in depth) is plotted along the CCD serial 
direction for different signal levels. The simulation includes an off -centred 
SBC and a saturation occurring at high signal level. 



where P refers to the parameter value in the buried channel and 
i^SBC to the corresponding value in the SBC. 

At high signal levels, saturation effects limit the linear growth 
of the charge density; while more charges can still be added, the 
maximum density cannot be overcome and the distribution expands 
in the three spatial directions. Saturation occurs when the charge 
density becomes larger than the doping concentration. Taking into 
account the presence of the anti-blooming drain we consider that 
locally, the charge density cannot exceed n sa t = ffsat/V L with 



V e = (27r) 3//2 a x a y a z . Following the method proposed by Linde- 
|gren| ( |199"8] ), the model for the saturation process is: 

_ n sat SV(x) 
neW "n sat + SV(x) 

where S' has to be adjusted to give the correct total charge S: 



(20) 



S 



with u = S'/S sa ,t and 5 sa t 
growth case: 



S — jjj n e (x,y, z) dx dy dz 

_ [J r°° _ 

V 7T J u 



(21) 



+ exp (r 2 /2) 
n sat x (2tt) 3/2 



(J x cr v <j z • In the linear 



(22) 



m< 1 and S ~ S sat u . 

In the saturation case: 

u > 1 and S ~ Ssat (4/3 v^F) (\nu) 3/2 . (23) 

The adjustment of S' is then analytically constructed to provide the 
right behaviour in the linear and saturation cases and a reasonable 
approximation for u w 1, i.e. the transition region: 

S' = S (1 + (S/S^ry 1 - 25 e (^s/M 2/3 . (24) 

Figure[3]shows an example of the evolution of the charge den- 
sity distribution as a function of the signal level as it is simulated 
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by the presented model. It illustrates the particular features that the 
model aims to reproduce. At low signal level the charge density dis- 
tribution shows a narrow profile in accordance with the expected 
SBC effect. While the signal level increases, the density profile 
widens and shifts to the centre of the pixel. It corresponds to the 
transition from the SBC regime to the BC regime in the case of a 
SBC implant not located at the centre of the pixel as depicted in 
Fig.[T]3 (which corresponds to the Gaia case). And finally one can 
clearly observe the saturation occurring for the large signal levels. 

2.4 CCD illumination history 

The CCD illumination history determines the CCD state (i.e. the 
trap occupancy level) prior to a star transit. Every star brighter than 
magnitude 20 will be observed ~ 80 times by Gaia (on average). 
During each observation of a given star, the satellite will be scan- 
ning the same part of the sky but from a different direction i.e. with 
a different orientation of the focal plane. Hence the occupancy of 
the traps in the CCD will be different prior to each transit of a 
given star. As a consequence the CTI effects (image location es- 
timation bias and charge loss) are likely to be different from one 
observation of a particular star to the next (even if the radiation 
damage to the CCDs were unchanged). Apart from the trapping 
process and the electron density distribution, the CCD illumination 
history is therefore another key element of the CTI effects mod- 
elling. The illumination history is determined by discrete events 
such as star transits and charge injections, which are directly re- 
produced during the simulation. Our simulations also account for 
a continuous optical background comprising light from unresolved 
stars and scattered light within the spacecraft. The light from the 
sky background constantly illuminates the Gaia CCDs at the level 
of rsj 42 e~ s _1 arcsec -2 (taking the contributions of both tele- 
scopes into account), i.e. each CCD pixel will on average receive 
~ 4 x 10 -4 e~ per pixel transfer step. In TDI mode the charges 
generated by the sky background are not only integrated but also 
continuously transferred from one pixel to another, to form a signal 
gradient along the parallel direction, the last light-sensitive pixels 
in the CCD effectively receiving about 2 e~. 

Brown & van Leeuwen (2009a) showed that such low levels of 
constant illumination noticeably modify the trap occupancy level. 
Indeed even if the local electron density distribution remains ap- 
proximately unchanged, the interaction time of these background 
electrons with the traps is effectively much greater than the one of 
a usual transiting source. To simulate the effect on the trap occu- 
pancy level, one can run the model for a certain amount of time 
with the background light as the only input signal. The star signal 
of interest is then inserted once an equilibrium trap occupancy level 
has been attained. However, our model is computationally very in- 
tensive and using it to simulate several minutes of CCD operation 
would take too much time. For this reason prior to each simulation, 
we determine the state of each trap using a random generator and 
considering the probability pmied oo, the probability for a particular 
trap to be filled by an electron from the background light, pmied oo 
is computed from eq. [12] assuming an infinite time of interaction: 

pmiedoo — pmied(t = oo) = — ^ — . (25) 

r r + r c 

In the simulations, the background contribution is also taken into 
account during the charge collection step by adding the background 
count rate to the expected number of collected photo-electrons be- 
fore invoking the Poisson random number generator (cf. Section 
|2.1| ). Note that the dark current can be taken into account in the 



same way. However in the Gaia case, the contribution of the dark 
current to the global background at operational temperature is not 
significant ( ~ 10 -7 e~per pixel transfer step), we thus ignore it in 
the following. 



3 RADIATION TESTS 

We describe the verification of our model against experimental 
data. These data were obtained by the industrial partners in the 
Gaia project. 

Several experimental studies were carried out on irradiated 
CCDs in order to evaluate the impact of CTI on Gaia's scientific re- 
quirements, to define the optimal operating temperature, to prepare 
the CCD calibration activities, and to elaborate a radiation dam- 
age mitigation strategy. Sira electro-optics acquired the first sets of 
CCD radiation test data (Hopkinson et al. 2005), focusing on the 
determination of trap parameters and CCD characterization. Later, 
Surrey Satellite Technology Limited (SSTL, formerly Sira) inves- 
tigated the potential difference between a CCD irradiated at room 
temperature and a CCD irradiated at 163K and kept at that temper- 
ature. SSTL concluded that the results obtained for CCDs irradi- 
ated at room temperature should be adequate for Gaia performance 
predictions within the usual experimental uncertainties ( Hopkinson 
[2008] >. 

Up to now the prime contractor for Gaia, EADS-Astrium, has 
performed three test campaigns on Gaia CCDs irradiated at room 
temperature with a radiation dose of 4 x 10 9 protons cm -2 (10 
MeV equivalent), and a fourth one is on-going. The experimental 
setup includes a translation stage which enables it to reproduce the 
star motion and hence to operate the CCD in TDI mode. The CCD 
was cooled and operated at constant temperature throughout the 
test campaigns. During the first campaign, Astrium concentrated 
on determining CTI effects on Gaia's astrometric measurements: 
the charge loss and the image profile distortion leading to a biased 
evaluation of the image location. The experiments were carried out 
as a function of stellar brightness (i.e. signal level) and background 
light level. The purpose of this first campaign was also to evaluate 
the viability of an artificial diffuse optical background source as a 
CTI mitigation device. The second campaign (RC2) allowed an al- 
ternative mitigation tool, Charge Injection (CI), to be thoroughly 
studied. Each of Gaia's CCDs contains a row of diodes and gate 
electrodes before the first pixel electrodes in the image section. A 
row of charges can thus be injected and transferred to fill the traps 
prior to a star transit. The CI level or the quantity of injected elec- 
trons is defined by the difference in voltage applied to the pixel and 
the gate electrodes. To study the influence of the CI parameters on 
the CTI effects, RC2 data has been acquired for different CI levels, 
durations (number of CI rows at a time), and delays (elapsed time 
between the CI and the first star transit), at different temperatures. 
The third campaign used a realistic sky-like illumination pattern to 
simulate the star transits instead of the uniform illumination grid 
from RC2. 

Not only astrometric tests were performed during these cam- 
paigns but photometric and spectrometric issues were also ad- 
dressed. The CCDs used by the red photometer and the radial ve- 
locity spectrometer (RVS) instruments on board Gaia differ slightly 
from the astrometric ones. They are based on the same architecture 
but are thicker devices and use a modified anti-reflection coating 
to enhance their quantum efficiency in the red wavelength band. A 
set of preliminary tests was performed during RC2 on such a de- 
vice irradiated with a lower radiation dose (2 x 10 9 protons cm -2 ), 
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followed by a more detailed study during the radiation campaign 3 
(RC3). This study included a mask mimicking realistic G2V stellar 
spectra and a very tight control over the background light. As a re- 
sult, the red-enhanced CCD was established to be more sensitive to 
radiation damage (most likely due to a greater depletion depth) and 
the expected shift and distortion of the spectral features induced by 
CTI were characterized as a function of stellar brightness and level 
of background light. 

The raw data acquired during these campaigns were made 
available to the Gaia Data Processing and Analysis Consortium in 
order to be re- analyzed and to support the CTI modelling efforts. 

To compare our model results with the experimental data, 
we relied mostly on the RC2 data (Georges 2008 ; Brown & van 
Leeuwen 2009b). We now describe in more detail RC2 experi- 
ments. The CCD irradiation scheme of the image section consists 
of three areas of similar size. The first area is non-irradiated, the 
second is irradiated at the level of 4 x 10 9 protons cm -2 (10 MeV 
equivalent) and the third, irradiated at a higher dose, was not used 
during testing. The serial register is non-irradiated. The reference 
and CTI affected data are obtained by translating the light source 
and a mask over respectively the first and second area of the CCD. 
The mask contains 50 x 22 (parallel x serial) holes, and the pro- 
jected stellar images are separated by 50 pixels in the parallel di- 
rection and 20 pixels in the serial direction. Each stellar image is 
binned along its serial dimension which necessitates some assump- 
tions in the input signal modelling (detailed in Section[4]). For each 
test at least 5 consecutive scans with identical configuration are per- 
formed. The mean time interval between the end of a scan and the 
beginning of a new one is 29.7 s. The first scan has a different illu- 
mination history from the others and is usually excluded from the 
analysis. We used two different astrometric tests out of the three 
that were performed during RC2. In the first test the CI delay and 
duration are kept fixed while the CI levels vary from ~ 4000 to 
^ 115 000 e~. Each injection level was tested at different temper- 
atures (from 163 K to 198 K) and for different illumination levels 
(corresponding to stellar magnitudes 13.3 and 15). In the second 
test the delay between the CI and the first star transit is varied from 
30 to at most 120 000 pixels (i.e. from ~ 29 ms to - 118 s) while 
the CI level is fixed to 20 000 e~ and the CI duration to 20 pixels. 
The test temperature was kept close to 163 K. The whole sequence 
of tests was repeated for different illumination levels (star magni- 
tudes 13.63, 15.29, 16.96, 18.65, 20.25). In the following section, 
we explain how we validated our MC model against these experi- 
mental tests. 



4 MODEL VERIFICATION AND COMPARISON TO 
EXPERIMENTAL DATA 

To validate the model we proceed in two steps. First the unit blocks 
of the model, the capture and release processes, are tested indi- 
vidually by comparing the results of Monte Carlo simulations to 
statistical analytical predictions derived from the equations in Sec- 
tion]^] The second step consists of a direct comparison between 
the model and a subset of the RC2 data. The variety of tests per- 
formed under different experimental conditions during RC2 allows 
us to separately verify the different features of our model. We first 
address the capability of our electron density distribution model to 
reproduce the CTI effects over a large range of signal levels and 
in particular to mimic the effect of the CCD supplementary buried 
channel. Subsequently we investigate how well the model repro- 
duces individual measurements acquired in TDI mode. 



4.1 Model comparison to analytical prediction 

During each transfer step of a simulation, it is possible to moni- 
tor the trap occupancy level # s im (or fraction of filled traps in the 
CCD): 



N 



(26) 



with iVmied, me number of filled traps and N, the total number of 
traps in the CCD. 

For a simple experimental configuration, one can analytically 
predict the evolution of the trap occupancy level # e x P as a function 
of time and compare it with # s i m to verify the reliability of the ba- 
sic steps in the Monte Carlo simulation. The capture module in the 
simulation is tested by simulating a CCD under a constant illumi- 
nation, which contains a unique trap species with a capture time 
constant r c and an infinite release time constant so that no electron 
release can occur. We also make all traps in the CCD empty at the 
beginning of the simulation, i.e. # s im(0) = 0. Under these condi- 
tions, the expected trap occupancy level # e x P is equivalent to the 
capture probability: 



#exp(t) = Pfilled : 



1 



(27) 



To compute r c analytically (eq. [8}, we assume a uniform electron 
density distribution, which we also implemented in our model for 
the purpose of this test. Similarly, to test the release module in the 
simulation we simulate a CCD in complete darkness, which con- 
tains a unique trap species with a release time constant, r r and an 
infinite capture time constant. All the traps are artificially filled at 
the beginning of the simulation so that # s i m (0) = 1. Under these 
conditions, one can write: 



#exp(t) — Pfilled — 1 ' 



(28) 



The number of traps in a CCD is finite. Thus, # s im and # e x P 
are necessarily different due to the discrete nature of the capture 
and release process (no fractions of electrons can be captured or 
released). Yet we want to assess how accurately our model repro- 
duces the expected value in these particular conditions. After a time 
t each trap in a CCD can only be empty or filled, with the proba- 
bility of occupancy given by pmied- Hence the simple simulations 
of trapping or release considered here constitute a Bernoulli trial 
where one counts the number of times n that a trap is filled or emp- 
tied after time t. The value of # s im is then given by n/N and will 
follow a binomial distribution with the following expectation value 
and variance: 

E(0 s im) — Pfilled — #exp , 
1 1 ( 29 ) 

Var(6>sim) = — pmied(l -pmied) = j^0 exp (l - 6> exp ) ■ 

We repeat the Monte Carlo simulation N run times in order to 
compute the mean trap occupancy level (0 S im) and its variance 
Var(# s i m ) and to verify whether these values are the same as the 
expected ones given in (eq.[29]). As can be seen from Fig. [4] for the 
capture and release tests the model reproduces the expected results. 
It is also interesting to note that even in this simplistic experimental 
configuration there is a clear difference (Fig. [4^) between the cap- 
ture profile obtained for a uniform electron density distribution and 
a more realistic distribution as described in Section [231 
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Figure 4. Tests of the capture (a) and release (b) processes. Top panels: comparison between the expected trap occupancy level (grey) and the mean trap 
occupancy level for 1000 Monte Carlo realizations for a CCD containing 100 traps (black). For the capture test (a), the dotted line depicts the results obtained 
for a uniform electron density distribution and the dashed line for a 3D Gaussian distribution including a SBC (eq.[T8). Bottom: the expected (grey line) and 
actual standard deviation (dotted black line, measured from the Monte Carlo experiments) for the trap occupancy levels, and the difference between them 
(black line). 



4.2 Model comparison to experimental data 

4.2.1 Fractional charge loss 

The charge loss induced by CTI is directly connected to the trap 
capture process, which is particularly sensitive to the electron den- 
sity distribution. As a consequence, reproducing charge loss mea- 
surements over a large signal range implies an accurate modelling 
of the electron density distribution. That is why we chose to extract 
fractional charge loss measurements from RC2 data to verify our 
particular approach to the electron density distribution modelling 
described in Section [231 

Thanks to its simple signal shape, it is particularly convenient 
to study the charge loss occurring in a charge injection (CI) line. 
An undamaged CI profile gives a constant signal with a mean value 
corresponding to the CI level (cf. Section [3}, whereas a damaged 
CI profile typically shows a strong electron deficit in the leading 
edge corresponding to electron captures and a slight electron sur- 
plus in the end of the profile as well as a trailing edge after the 
CI profile due to the release of electrons from the traps. To obtain 
the charge loss one needs to compare the damaged and undamaged 
profiles. RC2 data does not contain undamaged CI profile since the 
reference data were acquired without CI, it is therefore impossible 
to know with any accuracy, the number of electrons injected for 
a specific injection level. To compute the charge loss we then as- 
sume that the average number of electrons over the last 4 CI lines 
(CI duration = 20 pixels) would constitute an acceptable reference. 
Those last pixels undergo the least charge loss and therefore remain 
closest to the actual number of electrons injected per line. The sim- 
ulations, for which the true CI level is known, show that this as- 
sumption leads to a slight underestimation of the fractional charge 
loss at all signal levels. Hence, the last 4 CI lines are also used in the 
simulations to compute the CI reference level and subsequently the 
fractional charge loss to avoid any bias in our comparison. The frac- 
tional charge loss, formally the charge loss divided by the charge 
injection level times the CI duration, allows us to study the fraction 



of signal that is lost due to CTI in a CI profile as a function of the 
signal level. 

In Figure[5]the black crosses show the fractional charge loss as 
a function of the CI level extracted from RC2 data, we obtain results 
similar to Hopkinson et al. (2005) (for a direct comparison note that 
the devices were irradiated at different doses in the two studies). 
Fig.|5]reveals the complex structure of the pixel by showing a clear 
break (close to the nominal SBC full well capacity 1500 e~) in the 
increase of the fractional charge loss as the signal level diminishes. 

To reproduce these measurements with our model we simu- 
late a simplified version of the experimental setup. For each injec- 
tion level we simulate the operation of a CCD with a single pixel 
column of 4494 pixels, only 2 scans are performed (instead of 5), 
the CCD contains a unique trap species, no DOB is added, and the 
serial register is not simulated. The CCD is operated at the Gaia 
operational temperature (163 K), it is operated in TDI mode as in 
the RC2 tests. The selected electron density distribution includes 
the description of the SBC but no saturation. A charge injection is 
performed at the beginning of each scan. For each considered CI 
level, we perform this simulation with the same CCD (i.e. the same 
CCD characteristics and trap locations), we then measure the frac- 
tional charge loss following the same approach as in the test data 
analysis. The whole procedure is repeated 8 times and the charge 
loss measurements are then averaged so that the modeled fractional 
charge loss is independent from the noise due to the stochastic na- 
ture of the trapping processes. 

The model outcomes are then compared to an analytical fit to 
the experimental data and % 2 is used as a comparison criterion. 

2 = g (A (*Q- *(*,))' (30) 

2 = 

where A is the simulated fractional charge loss, $ the fitted frac- 
tional charge loss derived from RC2 data, Xi a particular signal 
level, and a the noise. For a particular signal level, the experimen- 
tal data is the mean of 4 measurements (cf. Section [3}. We assume 
that the standard deviation from this mean measurement (black bars 
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Figure 5. Comparison between the modeled fractional charge loss (grey) 
and experimental data extracted from RC2 (black) as a function of the signal 
level. 
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in Fig. [5} encompasses the experimental noise, the readout noise as 
well as the injection noise and is therefore equivalent to a. The set 
of parameters which minimizes our comparison criterion is found 
by using the downhill simplex minimization method (Nelder & 
Mead 1965 ). The free parameters in our simulation are the trap 
density, and the density distribution parameters: the SBC full well 
capacity, ?7bc, and ?7sbc- We introduce 77 in order to preserve in the 
electron density distribution the ratio between the signal confine- 
ment volume dimensions x macc . In this way the standard width of 
the distribution cr x in each direction corresponds to a fixed fraction 
of the predefined dimension of the signal confinement volume in 
that direction x macc : 

CT X = 77 X Xmax (31) 

To avoid local minima, we first randomly probe the parameter space 
and establish a x 2 map. The downhill simplex is then initialized 
with the set of parameters for which the x 2 is the smallest. The 
grey line and circles in Fig. [5] (left) give an example of the model 
results, representative of the quality of fit that can be achieved using 
this method. The simulation parameters for this particular fit are 
summarized in Table [4] The grey bars indicate 1-sigma deviation 
from the average over the 8 simulation iterations. Within the error 
bars, the model reproduces the experimental data over a wide range 
of signal levels (three orders of magnitude). The CTI mitigation 
effect of the SBC and the transition between the BC and SBC signal 
regimes are particularly well handled. At higher signal levels the 
performance of the model slightly degrades. We notice a deviation 
that would ultimately lead to underestimate the charge loss at very 
high signal levels. It is not possible to confirm this tendency as the 
charge injections with a signal level higher than 17 ke~ could not 
be studied due to saturation of the output amplifier. 

A comparison between our analytical description of the elec- 
tron density distribution and the more detailed model of a Gaia 
pixel (Seabroke et al. 2009 2010 ) shows that in the serial direc- 
tion the Gaussian approximation may be too crude at high signal 
level, where the distribution saturates quickly and takes a box-like 
shape. However, when using Seabrokes model we were not able 
to reproduce satisfactorily the fractional charge loss measurement. 
The simulations deviated significantly from the experimental data 
in particular during the transition between the BC and SBC sig- 



Table 4. Simulation parameters and goodness-of-fit for the model result 
example. The fitted parameters are indicated with a 



nal regimes. One possible explanation is that the single e2v Gaia 
CCD in question has a non-nominal pixel architecture due to man- 
ufacturing tolerances (Seabroke's model uses nominal e2v design 
values). 



4.2.2 Astrometric images 

The verification of the elementary units of the model being satis- 
factory as well as the validation of our electron density distribution 
approach, we are now interested in estimating the capabilities of 
the model to reproduce the CTI induced distortion of stellar images 
acquired in TDI mode. 

Thanks to the second test of RC2 (cf. Section [3}, the model 
performance can be investigated for different stellar brightnesses. 
To do so we compare the model outcomes with the stellar images 
acquired in the irradiated part of the CCD. This comparison is per- 
formed at the sub-pixel level thanks to damaged over- sampled stel- 
lar images built using the multiple scans performed for each set of 
experimental parameters (temperature, star brightness, CI delay). 
The model input signal must be representative of the CCD illu- 
mination conditions. Therefore, we used the accumulated data in 
the non-irradiated part of the CCD to create a reference undam- 
aged image profile. In order to stay as close as possible to the orig- 
inal experimental set-up, the transit of the investigated star is sim- 
ulated in two dimensions, i.e. over a virtual CCD with 12 pixel 
columns of 4494 pixels. This requires a two-dimensional input. 
However during RC2, the images were binned in the serial direc- 
tion to recreate in-flight conditions such that the reference curves 
are LSFs (Line Spread Function, the PSF integrated in the serial 
direction). Since we have no information on the two-dimensional 
PSF, to generate a PSF from the original reference curve we as- 
sume that the profiles in the parallel and serial directions are the 
same: P (x, y) = L (x) x L (y), where L is the undamaged refer- 
ence curve, x and y the positions in pixel respectively in the parallel 
and serial directions, and P the model input image. The integrated 
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flux of the reference curve is scaled to produce an input image for 
each level of illumination (or artificial star magnitude). 

To perform a direct comparison between the model and 
the over-sampled damaged profiles, we first extract the sampling 
scheme specific to each of the damaged profiles. We then apply it 
to the PSF generated from the reference curves to create the model 
input signal. In this way the required number of simulations for a 
single set of parameters to generate an n times over- sampled sim- 
ulated damaged profile is n. Once the simulations are completed, 
the individual predictions are binned in the serial direction and each 
data point is then placed at the correct sub-pixel position according 
to the original sampling scheme so as to form an oversampled sim- 
ulated damaged profile. 

We present the results of this detailed comparison (Fig. [6] Ta- 
ble|5]> for two different cases: (i) the model is first used to reproduce 
a single image profile (i.e. for a unique set of experimental param- 
eters: temperature, CI delay, and magnitude); (ii) and then, with a 
single set of simulation parameters, the model is used to simulta- 
neously reproduce a set of damaged image profiles with different 
magnitudes. 

When fitting to an individual profile (Fig.[6](a)), % 2 is used as 
our goodness-of-fit criterion: 



X 



s-i 

E 



(X(x l )-N(x l )Y 



(32) 



where A is the simulated damaged profile, N the RC2 profile, x% a 
particular sub-pixel position, S the total number of data points, a, 
the noise, is considered to be equivalent to the quadratic sum of the 
photon-noise and the readout noise r. The photon-noise is assumed 
to be Gaussian with a standard deviation of y/N (an approximation 
of the Poisson statistics for large photon counts) and r is assumed 
to have the constant value of 4.8 electrons: 

a 2 — N ( Xi ) + r 2 (33) 

When fitting to a set of damaged profiles (Fig.[6](b)), \ 2 is altered 
in order to avoid any fitting bias towards the brightest magnitudes 
and the most over-sampled profiles. The new comparison criterion 
g is thus defined as follows: 
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(34) 



where F is the total integrated flux. 

In each case the free parameters of the simulation are the trap 
parameters (p the density, a the capture cross-section, and r the re- 
lease time constant). The number of trap species is fixed to 2. The 
electron density distribution parameters are fixed to the values sum- 
marized in Table [4] Compared to the fitting of the fractional charge 
loss measurement, fitting astrometric measurements is a complex 
task: 

(i) The parameter space is larger. 

(ii) The parameter space is more degenerate; a single image pro- 
file does not set high constraints on the trap parameters. 

(iii) Both the capture and release processes are involved. 

(iv) The illumination history plays an important role. 

(v) As we will explain in the following, the number of experi- 
mental uncertainties is higher as well as the number of assumptions 
required by the data processing and the simulation. 

(vi) The fitting also necessitates a high number of two- 
dimensional simulations of stellar transits, which is computation- 
ally very intensive, and sets a high requirement on the minimization 
procedure efficiency. 
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Table 5. MC model fitting parameters corresponding to the examples shown 
in Figs, ©a) and b). Note that X 2 ed , the reduced x 2 , indicates here the 
goodness-of-fit. In the three presented cases, a short and a long release time 
constant species are needed to explain the experimental data. 



The preliminary % 2 maps also appeared to be hard to characterize 
as they contained a lot of local minima, none of which presented 
a satisfactory agreement with the data. In order to better sample 
the parameter space we decided to perform the first step of our 
minimization procedure by means of an evolutionary algorithm^] 
We applied two evolutionary mechanisms, mutation and cross-over, 
with optimal occurrence probabilities on a restricted initial popula- 
tion of parameter sets and let it evolve towards smaller comparison 
criteria. After a limited amount of generations, we initialized the 
downhill simplex algorithm with the set of parameters for which 
our comparison criterion was the smallest. As previously men- 
tioned, simulating two-dimensional stellar transits is a relatively 
time-consuming process. For long charge injection delays (> 20s) 
one can neglect the effect of the charge injection release trail, the 
star then crosses a CCD with almost all the traps empty (from short 
to fairly long release time-constants). This enables us to drastically 
reduce the number of simulated TDI steps. Additionally the stellar 
transits acquired in these conditions offer the advantage of present- 
ing important CTI effects that set a higher constraint on the fitted 
parameters. We thus selected the subset of RC2 damaged profiles 
with a CI delay of ~ 30 s to perform our comparison. Finally, al- 
though the profiles were ultimately compared at the sub-pixel level, 
in order to decrease the total number of simulations, only one sam- 
ple per pixel was simulated during the fitting procedure Figure [6] 
(a) presents an example of a simulated damaged profile (circles), 
this profile is representative of the best-fitting achievable in these 
particular conditions at a magnitude of 13.67 and should be com- 
pared to the black line profile (RC2 data). As can be seen from the 
residuals normalized by the photon-noise, the simulation is in fairly 
good agreement with the data over the whole profile. However we 
note a slightly larger disagreement in the image leading edge. This 
is not surprising as for bright images the wings still contain a fairly 
large amount of electrons (e.g., 800 e~in the first simulated sam- 
ple of the leading edge), which may play an important role in the 
self-illumination history by filling a number of traps and mitigat- 
ing the CTI effects in the image core. Hence the fit may be further 
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Figure 6. (a) Top left: Comparison between the RC2 test data (black line) as reduced by Brown & van Leeuwen ( 2009b} and the model simulation (black 
circles) for a particular signal level and CI delay (Magnitude: 13.63; CI delay: ~ 30 s). The charge transfer direction is from right to left. The amplitude of the 
radiation damage can be appreciated by comparing the damaged profiles (black line and circles) to the input signal (dotted line) based on the non-irradiated 
test data. Bottom left: difference between the simulation and the experimental data normalized by the photon-noise. The quality of fit is similar in the profile 
core and wings. Note the clustering of the sampling: each data point within a pixel corresponds to a different scan. 

(b) Top right: Comparison between the RC2 test data (crosses) and the model simulations (circles) for several signal levels and a particular CI delay (Magnitude: 
13.63 (black), 15.17 (dark grey), 17.03 (light grey); CI delay: ~ 30 s). The model simulations are obtained for a single set of parameters (Table|5j, which has 
been optimized to find a reasonable agreement with the test data. Note that this time the ordinate scale is logarithmic in order to facilitate the visual evaluation 
of the quality of fit for each magnitude and in the core as well as in the wings of the image profiles. The charge transfer direction is from right to left. Bottom 
right: Comparison between the resulting experimental (continuous line) and simulated (dashed line) charge loss as a function of magnitude. The dotted line 
corresponds to the resulting charge loss for the simulation parameter obtained when fitting to magnitude 13.63 only. 
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Figure 7. Top: Comparison between the RC3 test data (black line) as reduced by Hall ( 2010 1 and the model simulation (grey line) for a star of magnitude 12.5 
(No CI, and lowest level of background light (< 1 e~ /pixel)). The charge transfer direction is from right to left. The reference averaged spectrum (acquired 
in the non-irradiated part of the CCD) was used as input signal (dotted). Bottom: The resulting charge loss normalized by the photon-noise for the simulation 
(grey) and the test data (black). 
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improved by simulating a profile wider than 15 pixels. The best in- 
dividual fits obtained at fainter magnitudes (cf. example in Table[5]) 
shows a slightly better overall agreement. The relative amplitude of 
the CTI effects in the core compare to the wings is very sensitive to 
the input image shape in the CCD serial direction. The assumption 
made in order to build a two-dimensional input signal may thus also 
limit the ultimate goodness-of-fit. 

Figure [6] (b) shows the resulting fits at three different magni- 
tudes (13.67, 15.17, 17.03) obtained for a single set of simulation 
parameters (cf. Table [5}. The overall agreement is remarkable. As 
can be seen in the bottom part of the plot, the amplitude of the CTI 
effects is qualitatively well reproduced as a function of the signal 
brightness. And as expected the charge loss peak slowly shifts to- 
wards the image core with fainter signal. However one should no- 
tice that when a single set of trap parameters is fitted to several pro- 
files at different signal levels the resulting goodness-of-fit for each 
individual profile decreases, see Table [5] This is also illustrated in 
the bottom part of Fig.|6](b), where the resulting charge loss curves 
at the magnitude 13.67 for the individual fit case (black dotted 
line) can be compared to the multiple magnitudes fit case (black 
dashed line). This difference in goodness-of-fit may be ascribed to 
variations in the experimental conditions between two tests with 
different signal level (slightly different temperatures and illumina- 
tion background). Also the intervals between the tests were several 
days. In between two tests the CCD was stored at room tempera- 
ture. Astrium later acknowledged that this may have resulted in a 
change of the CCD state, as discrepancies in test results were ob- 
served for similar experimental conditions. As already discussed 
in[Prod'homme et al. (2010), this may set a limit on the ultimate 
goodness-of-fit achievable with a single set of trap parameters to 
several profiles by any model. As a final remark on this part of 
the comparison, we would like to state that even if the trap species 
parameters found as a result of the fitting procedure could be asso- 
ciated with known trap species from the literature, the uncertainties 
in the experiments are too large and the assumptions in the sim- 
ulation process too numerous to conclude that such species were 
indeed present in the tested CCD. If one wants to infer trap species 
parameters by using this model to reproduce test data, the experi- 
ment should be carefully designed to serve this purpose. Artificial 
charge injections would be particularly suited for such an experi- 
ment as one can infer trap densities and capture cross-sections from 
the charge loss that occurs in their profile and release time constants 
from their trailing edge. We would in particular recommend: (i) to 
repeat the tests at different temperatures in order to break the po- 
tential degeneracies in the time release constant space; (ii) to use 
different levels of charge injection from a few electrons to the pixel 
full well capacity, in order to make sure that no trap is missed as 
well as to set the highest constraint possible on the charge density 
distribution model over a complete signal range. 

4.2.3 Radial velocity spectrometer images 

The radial velocity of stars will be measured on-board Gaia thanks 
to a medium resolution spectrometer ( |Katz et aL]|2004| ) that will 
enable the analysis of the Doppler shift of the stellar spectral lines, 
in particular the Ca II triplet. The spectrometer is composed of a 
series of prisms and a diffraction grating that will disperse the stel- 
lar light before it reaches the red-enhanced Gaia CCDs. The com- 
bined effect of TDI and high light dispersion will result in a very 
faint spectral signal, down to a few electrons per spectral feature. 
The signal-to-noise ratio for a single measurement is expected to be 
very low (e.g., 7 for a star of magnitude 16), and, at the faint end, 



the radial velocity measurement will rely on the co- addition of mul- 
tiple spectra to recover the spectral features. Until recently CTI had 
never been studied at such low signal levels, and it was uncertain if 
one could indeed recover spectral features by co-adding damaged 
spectra. RC3 has been designed to address this particular issue by 
the use of a mask manufactured to reproduce in detail the spectrum 
of a G2V star. To illustrate the capacity of our model to reproduce 
the CTI effects in TDI mode at extremely low signal levels (< 20 
e~) we use the RC3 test data and compare the model outcomes to 
the stellar spectra acquired in the irradiated part of the CCD. 

The model input signal is built from the accumulated data in 
the non-irradiated part of the CCD and the transits are simulated 
over a virtual CCD containing 6 pixel columns of 4494 pixels. Dur- 
ing the tests the spectra were binned over 12 pixels in the serial 
direction, we thus again had to assume the serial signal profile to 
perform two-dimensional simulations. We assumed the same par- 
allel spectral profile in every pixel column but with different in- 
tensities. The fractions of the total intensity for each pixel column 
were chosen according to the serial profile of the simulated spec- 
trum provided by the Gaia Instrument and Basic Image Simulator, 
GIBIS (0.04, 0.14, 0.26, 0.29, 0.20, 0.07). This time we do not per- 
form any optimization of the trap parameters but we selected two 
trap species from |Seabroke et al. (2008) (divacancy and unknown: 
pi, 2 =1 trap/pixel, cri j2 = 5 10~ m , n — 20 s, and T2 = 80 
ms). Once again we fixed the electron density distribution parame- 
ters to the values summarized in Table [4] 

In these conditions, one can obtain a spectrum such as the one 
presented in Fig.[7]for a star of magnitude 12.5 (grey line). In a sim- 
ilar fashion to the depicted experimental damage spectrum (black 
line), 15 transits were averaged to obtain the simulated spectrum. In 
both the RC3 data and the model outcome, one can notice that even 
at such low signal level the CTI effects do not completely wash 
away the spectral features; they can be recovered by co-adding 
several spectra. However, as can be seen in more detail from the 
blown-up portions of the presented spectra, CTI lowers the contin- 
uum and affects the absorption features by increasing their width 
and shifting them. These effects lead to a significant decrease of 
the overall spectral contrast. In the bottom part of the figure, we 
show the charge loss normalized by the photon-noise. The model 
reproduces remarkably well the first significant charge loss bump 
(only few electrons) occurring at the leading edge of the spectrum 
(left). The electron release that occurs at the trailing edge after the 
spectrum transit (right) is also well reproduced. The simulated av- 
eraged profile is binned in the serial direction over 6 pixels instead 
of 12 during the tests, as a consequence the simulation is slightly 
noisier than the RC3 data. It is thus hard to distinguish, but within 
the spectrum itself one can notice that the model does widen and 
shift the spectral features (see the blow-ups). The amplitude of the 
shift and increase in width is not quite matched though, this could 
be imputed to a difference in the release time constant for the trap 
species. 



5 CONCLUSION 

We have described a physical Monte Carlo model that simulates 
CTI effects induced by radiation damage in astronomical CCDs at 
the electrode level. This model has been elaborated in the chal- 
lenging context of ESA's Gaia mission. The operating conditions, 
the extreme astrometric requirements and a novel CCD pixel ar- 
chitecture, necessitated the development of a new approach in the 
representation of the electron density distribution as well as a more 
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detailed description of the trapping probabilities. These new fea- 
tures have been combined in a comprehensive simulation of CCD 
charge collection and transfer at the electrode level. In order to ver- 
ify the model predictions, we first validated the unit blocks of the 
simulation and then proceeded to a detailed comparison between 
the model and experimental test data. We showed that the model 
is able to accurately reproduce the CTI effects for a wide range of 
signal levels down to a few electrons, hence validating our electron 
density distribution representation. We finally validated the global 
model operation by assessing the model capability to reproduce the 
CTI induced distortion on stellar images and spectra acquired in 
TDI mode at different magnitudes. 

The model elaboration contributed greatly to our present un- 
derstanding of CTI effects. In particular, by showing that the im- 
plementation of a density driven approach of the electron packet 
growth enables the reproduction of experimental data and that to be 
successful in modelling the CTI effects at very faint signal levels, 
no detail should be neglected. The simulations should be as real- 
istic as possible, down to the transfer of electrons at the electrode 
level and the simulation of each individual trap. 

Due to the complexity of the CTI effects (including the vary- 
ing illumination history) and the extreme accuracy required, in par- 
ticular on the estimation of the stellar image location, one cannot 
apply a conventional correction of the raw data. Even if this were 
possible, the process would likely be unstable and lead to poorly 
understood error propagation. Instead of a direct correction, the 
Gaia Data Processing and Analysis Consortium (DPAC) is devel- 
oping a scheme which relies on a forward modelling approach that 
enables the estimation of the true image parameters from a dam- 
aged observation (e.g., Prod'homme 2010). In this scheme each 
observation is ultimately compared to a modelled charge profile in 
which the distortion of the CTI-free image (PSF) will be simulated 
through an analytical CTI model. Currently the Monte Carlo model 
described in this paper is used to generate large synthetic data sets 
of both damaged and damage-free observations. The simulations 
are used to re-assess the final performance of the mission, taking 
the effects of radiation damage into account, as well as to verify 
the DPAC CTI mitigation scheme ( |Prod'homme et al.|201l] [HoITl 
|etal.|2011) . 

Although developed in the Gaia context, we tried to keep the 
model as general and as flexible as possible. It can be used to sim- 
ulate any kind of measurements performed with a CCD operated in 
imaging mode or TDI mode. Different clocking schemes can be ap- 
plied and the description of the electron density distribution should 
be flexible enough to simulate different pixel architectures. The use 
of this model is particularly relevant in the frame of space experi- 
ments that aim at very accurate measurements at low signal levels. 
ESA's Euclid mission, the astrometric measurements performed on 
board HST or future X-ray missions sent to L2 may benefit from 
the use of such a model to evaluate the impact of radiation damage 
on their performance budgets. 

A running version of the model as well as a brief documenta- 
tion and a few examples are readily available at jhttp : / /www . | 
|strw . leidenuniv . nl/ ~prodhomme7 cemga . php as part 
of the CEMGA java package (CTI Effects Models for Gaia) de- 
veloped at Leiden Observatory. Please contact the authors for more 
information on how to use the package. 
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